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We present simulation results addressing the dynamics of a colloidal system with attractive in- 
teractions close to gelation. Our interaction also has a soft, long range repulsive barrier which 
suppresses liquid-gas type phase separation at long wavelengths. The new results presented here 
lend further weight to an intriguing picture emerging from our previous simulation work on the 
same system. Whereas mode coupling theory (MCT) offers quantitatively good results for the de- 
cay of correlators, closer inspection of the dynamics reveals a bimodal population of fast and slow 
particles with a very long exchange timescale. This population split represents a particular form 
of dynamic heterogeneity (DH). Although DH is usually associated with activated hopping and/or 
facilitated dynamics in glasses, the form of DH observed here may be more collective in character 
and associated with static (i.e., structural) heterogeneity. 
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I. INTRODUCTION 

Colloidal systems have an important experimental role 
as model materials for studies of the glass transition 1 - 2 . 
In particular, it is possible to approach the limit of hard 
sphere particles, and then deviate from this in a con- 
trolled manner by the addition of interactions, such as the 
depletion interaction mediated by added polymer 3 . Re- 
cently this has allowed the role of short range attractions 
to be probed experimentally^; that work has confirmed 
a scenario of re-entrant melting (on addition of a weak 
short ranged attraction) followed by re-freezing into an 
attraction-driven glass. This scenario was first suggested 
by MCT (mode coupling theory) calculations^^. Its 
confirmation in experiments represents strong evidence 
that MCT, while admittedly deficient in its neglect of 
activated local dynamics2*i2*ii, does capture some of the 
important collective dynamics of the colloidal arrest tran- 
sition. 

Within MCT, these collective dynamics are dominated 
by colloidal interaction forces (repulsive caging or at- 
tractive bonding) rather than details of the local diffu- 
sivity. For this reason, when simulating such systems, 
many-body hydrodynamic coupling is frequently ignored, 
replacing the full hydrodynamic interactions with local 
drag (Brownian dynamics: BD). Moreover, if the local 
dynamics is truly unimportant, one may even replace 
Brownian dynamics with Newtonian (i.e., molecular dy- 
namics: MD) which further saves computer time. In ef- 
fect the solvent has then been replaced by a vacuum. This 
brings the simulated system very close to that of atomic 
or molecular glasses, although the interaction parame- 
ters are typically somewhat different: a shorter range of 
attraction is common in colloids. 

Here we report various MD simulation results on model 
colloidal systems with short-range attractions, close to 



the arrest transition (gelation). These extend the work 
presented by ourselves in a series of earlier papersiS^i 3 ^ 
which revealed a somewhat paradoxical picture of dy- 
namics in these systems. We found 12,13 , that MCT is 
quantitatively successful at predicting correlation func- 
tions and other ensemble-averaged properties close to ar- 
rest; this finding offers strong evidence for the physical 
assumptions of collective rather than activated dynam- 
ics underlying MCT. Yet, when we looked at individual 
particle motions^, the behavior is far more complicated 
than we expected, revealing a strongly bimodal distri- 
bution of particle displacements. These can be resolved, 
over prolonged time scales, into two separate populations 
of fast and slow particles (with a very sluggish exchange 
between these). Moreover the fast and slow particles 
are highly correlated in space; at any instant, the dy- 
namics resembles a near-frozen gel of slow particles with 
fast particles moving around and through channels in the 
gelii Part of this picture may result from a long-range 
repulsive contribution to the chosen interaction potential 
whose role is explained below; nonetheless, the results 
can shed important light on the nature of collective and 
local motions in arresting systems. 

The behavior that we have observed represents a spe- 
cific form of dynamical heterogeneity (DH), related to, 
but distinct from, the DH found for glasses arrested by 
cagingi^iS*iLi2*iS. The latter observations of DH tend to 
support a 'facilitated dynamics' picture of the glass tran- 
sition in which regions of high mobility interact nonlin- 
early to give a super-activated temperature dependence 
of the overall mobility. Models of this type have been 
studied for several years 2 ^ and lately raised to a higher 
level of predictive power by Chandler and coworkers^!. 
However, these approaches generally model a dynamics 
in which a frozen matrix (often represented by a lattice) 
is relaxed progressively as the heterogeneities visit differ- 
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ent neighborhoods. Such models do not yet capture the 
rather rich picture observed in our simulations; for exam- 
ple, fast particles are found preferentially at the surface 
of the frozen component and individual bonds between 
fast particles have shorter lifetimes than bonds between 
slow ones^. These features probably cannot be under- 
stood within a lattice DH model but would require ex- 
plicit treatment of off-lattice dynamics in the presence of 
short range bonding. 

An important question is why the presence of such het- 
erogeneities does not destroy the agreement with MCT 
found in the ensemble-averaged dynamics, as measured 
in correlation functions^ 3 -. We do not yet have a clear an- 
swer to this question. However, we note that the hetero- 
geneous dynamics observed, once the fast and slow com- 
ponents are identified, is correlated with structural het- 
erogeneity; and this is (at least partly) collective in origin. 
The heterogeneity is therefore not simply attributable to 
the activated 'local hopping' channel which was long ago 
recognized to be missing from MCT2*i£. Indeed, some as- 
pects of this population- level dynamic heterogeneity may 
have little to do with the physics of activation (although 
that does seem to control a slow exchange between pop- 
ulations that we observe; see Section lill Bl below') . 

If the observed dynamic heterogeneity reflects in part 
collective rather than activated processes, then its pres- 
ence in a system whose overall properties (as measured 
by decay of globally averaged correlators) are well de- 
scribed by MCT is more easily explicable. Such dynamics 
could arise, for example, by an interplay between MCT- 
like arrest and incipient liquid-gas like phase separation. 
Macroscopic phase separation is prevented, in our simula- 
tions, by a barrier in the potential^. In real life, even for 
systems quenched beyond the relevant spinodal, it could 
be prevented instead by the extremely hindered sepa- 
ration dynamics that arise if one of the two coexisting 
phases is nonergodic. Incipient fluctuations of the same 
physical character should also exist close to the spinodal 
but within the stable phase. (In our own work we do 
not see microphase separation but do not rule out an 
incipient or frustrated form of this.) 

The interplay of phase separation and arrest has re- 
cently been discussed for colloidal systems at relatively 
low density, where some sort of iterative renormalisation 
of the MCT parameters is required^ 2 .. At higher colloid 
density, as we study here, MCT could in principle cap- 
ture some of this interplay unaided; after all, the input for 
MCT is the static structure factor S(q) which contains, 
as q — > 0, information about proximity to the liquid-gas 
spinodal. (In our system where a barrier in the potential 
suppresses this spinodal, S(q) contains instead informa- 
tion about any incipient microphase separation through 
the 'pre-peak' at low q.) However, we emphasise that 
in other aspects of the observed dynamics, particularly 
connected with interchange between fast and slow pop- 
ulations, MCT could well offer only limited insight, as 
discussed below. A somewhat complicated picture in- 
volving activation and/or facilitated dynamics, perhaps 



coupled to structural inhomogeneity and collective dy- 
namical modes, may ultimately be required to model the 
system studied here. 

In what follows, we first recall in Section [H] the details 
of our simulations. We then present new data that offers 
further quantitative support for the picture of dynam- 
ics close to colloidal gelation that was outlined above. 
In Section fill Al we add to our MCT-inspired analysis^ 
by presenting new results for the coherent autocorrela- 
tor, its wavevector-dependent nonergodicity parameter 
f q , and its wavevector-dependent terminal time T q . In 
Section lill Bl wc add to our population-based analysis of 
DHi^ with additional quantitative measures of the differ- 
ing local environments of fast and slow particles, based 
on the number of neighbors of each type and on measures 
of percolation. In Section Hvl we give our conclusions. 



II. SIMULATION DETAILS 

We have performed computer simulations of a sys- 
tem composed of soft core polydisperse particles, with 
a short-range attraction given by the Asakura-Oosawa 
potential 3 , modeling a mixture of colloids with non- 
adsorbing polymers. Phase separations have been 
avoided, in order to have full access to the whole composi- 
tion space. Crystallization is inhibited by polydispersity, 
and liquid-gas demixing by a long range repulsive bar- 
rier. The gel transition is generally approached from the 
fluid side, by increasing the polymer fraction within the 
fluid phase. 

The total interaction has three parts, a core repulsion, 
a short range attraction and a repulsive barrier. The core 
repulsion is given by: 



V sc (r) = k B T ( — ) 
\a 12 J 
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where an = (a± +02), with a\ and a-i the radii of the in- 
teracting particles. Particle sizes are distributed accord- 
ing to a flat distribution of half-width S = 0.1a, where a 
is the mean radius; a^ 6 [0.9a, 1.1a]. The attraction in- 
duced by the polymers, extended to take polydispersity 
into account, reads&2&24: 



VaoM = -k B T4> P 



(^l) 3 -|(^l) 2 + 7^ 
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for 2(a i2 + £/5) < r < 2(a X2 + £) and V A a = for larger 
distances. Here, 77, = fj = (771 + 772 ) / 2 , and cj> p is 

the volume fraction of the polymer. The range of the po- 
tential is given by £, the polymer size, and its strength is 
proportional to <p p . At short distances, for computational 
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efficiency this potential is replaced by a steep parabola 
with the minimum at r = 2a 12 ; this ensures that the 
total interaction potential retains a quadratic minimum 
very close to r = 1a\2. 

The long range repulsive barrier inhibits liquid-gas sep- 
aration by raising the energy of the dense phase. It has 
the form: 



Vbar(r) — ksT ■ 



r — r\ 
So - n 



r — r± 

ro - n 



(3) 



for ro < r < n and Vb ar (r) = otherwise. The limits 
of the barrier were set to ro = 2(ai2 + £), and r\ = 4a, 
and its height is only lfcsT. This energy is equal to the 
attraction strength at 4> p = 0.0625. As described in the 
introduction this leads to a low-g peak in the structure 
factor but this is relatively mild (always lower than the 
main peak) and various tests for microphase separation 
gave negative results 1 ^. However, as previously discussed 
we do not rule out some effect of frustrated or incipient 
microphase separation on the dynamics. 

The resulting total interaction potential, Vtot = V sc + 
Vao + Vbar , is analytic for all distances and can be used 
with simple molecular dynamics algorithms. Polydis- 
persity in size causes differences in the interaction be- 
tween different pairs of particles; apart from the obvi- 
ous differences in the core radii the attractive well depth 
also increases with particle size, the maximal variation 
beingi&ii 1.25 k B T. 

The simulated system is composed of 1000 particles, 
length is measured in units of the average radius, o, and 
time in units of y/4a 2 /3v 2 , where the thermal velocity, 
v, was set to y4/3. Equations of motion were integrated 
using the velocity- Verlet algorithm, in the canonical en- 
semble (constant NTV), to mimic the colloidal dynamics, 
with a time step equal to 0.0025. Every n t time steps, the 
velocity of the particles was re-scaled to assure constant 
temperature. No effect of n* was observed for well equili- 
brated samples. Equilibration of the samples was tested 
by the trends of the energy of the system and other or- 
der parameters 1 ^, and by the time-translation invariance 
of the correlation functions. The range of the attrac- 
tion is set to 2£ = 0.2a, and the control parameters are 
(i) the density of colloids, reported as volume fraction, 

(f) c = |7ra 3 ^1 + (|) 2 ^ n c , with n c the colloid number 

density, and (ii) the polymer volume fraction, <p p , which 
measures the attraction strength. 

Because of the dynamical heterogeneities observed, 
care is needed in gaining adequate statistics over quasi- 
independent realisations. For every state point, differ- 
ent systems were equilibrated by quenching from hard 
spheres to the desired <f> p . For the slowest evolving states 
(cf)p = 0.405, 0.41, 0.415 and 0.42) ten different quenches 
were performed and 100 different time origins were con- 
sidered for computing correlation functions, all of which 
lay beyond the initial equilibration period. The reported 
correlation functions for these states are thus, in effect, 




FIG. 1: Snapshot of the system close to the gel transition: 
(j> c = 0.40 and <j> p = 0.42. 



averaged over 1000 independent realizations. The equili- 
bration period was as long as 5 x 10 4 time units in the 
slowest cases; this was enough to give time-translation in- 
variant correlators but nonetheless may not be long com- 
pared to the timescale of exchange of particles between 
fast and slow populations 14 . The issue of whether sam- 
ples are "fully" equilibrated is therefore partly a matter of 
interpretation - as is often the case in systems approach- 
ing an arrest transition, and always the case within the 
gel phase itself. 



III. RESULTS AND DISCUSSION 

The system under study has been specially devised to 
allow us full access to the whole (<p c ,<p p ) plane, eliminat- 
ing both crystallization and liquid-gas demixing. How- 
ever, the interplay between the long range repulsive bar- 
rier and the strong short range attractions induces a 
somewhat heterogeneous structure in the system, with 
voids and tunnels (see Fig. QJ, which minimizes the po- 
tential energy. Such a heterogeneous system is charac- 
terized by a peak at low (but finite) wavevector in the 
structure factor, S(q). In our system, such a pre-peak 
forms as the attraction strength is increased, as observed 
in Fig. El showing the build up of the heterogeneous 
structure. 

A low angle peak has been also noticed in experiments 
on colloidal gelation, moving to lower wavevectors as time 
proceeds, until it finally arrests 2 ^^. Despite the appar- 
ent similarity between both peaks, their origin is quite 
different; whereas the experimental one is due to the 
arrest of liquid gas demixing by gelation, the peak in 
the simulations is an equilibrium feature, caused by the 
long range repulsion. Nonetheless, the effect of these two 
mechanisms on dynamics within the final structure could 
be similar, and we believe that MD studies of the sys- 
tem with barrier offer the closest numerical analogue cur- 
rently available for such experiments. (If no barrier were 




FIG. 2: Structure factor for different states on the isochore 
4> c — 0.40. Approaching the gel transition: 4> p = 0, 0.20, 0.35 
and 0.42. 

used, because the simulation is small compared to a real 
system, substantial phase separation on the length scale 
of the simulation box would occur quite rapidl y 13 ! 27 .) 

Under certain circumstances, long range repulsive bar- 
riers are known to cause microphase separation-— -—i-l. 
In that case, the regions of high and low particle density 
arrange in a periodically ordered pattern, causing one or 
more Bragg peaks at low q in equilibrium. In our case, 
the barrier is not strong enough to cause this microphase 
separation, as observed by the disordered structure of the 
voids and tunnels, and the height of the prepeak. The 
system is thus an amorphous fluid, with frustrated liquid- 
gas separation, but will possibly separate in microphases 
at higher attraction strength. This is similar to other dis- 
ordered phases showing incipient microphase separation, 
such as microemulsions 31 . 

In the following sections we make use of the density au- 
tocorrelation function, or normalized intermediate scat- 
tering function: 



^ = n \ £ £ 6XP {tq ' M<) ~ Tjim ) /S(q) 

\ 1=1.7 = 1 / 

(4) 

where N is the number of particles in the system, and 
q is the wavevector. The self part of it, & q {t), is cal- 
culated by restricting the double summation to the case 
i = j. The brackets imply ensemble average, which we 
have calculated using different time origins as well as dif- 
ferent samples, as explained above. 



A. MCT analysis 

Within mode coupling theory, the main region of 
the arrest line that represents an attraction-driven 
(gelation-like) transition corresponds to a so-called A2 
bifurcation*. The same is true for a standard, caging 
driven transition; but some differences with respect to 




FIG. 3: Time-rescaled density autocorrelation functions for 
(f> p = 0.39, 0.40, 0.405, 0.41, 0.415 and 0.42. 

the latter are predicted in the non-universal parameters 
of the bifurcation, such as the nonergodicity parameter. 
Such differences arise from the different driving mecha- 
nisms: caging by steric hindrance for the repulsive glass, 
and formation of a bonded network in attractive glass. At 
high density and attraction strength, where the repulsive 
and attractive glass lines meet, a high order singularity is 
found for short range potentials, characterized by a log- 
arithmic decay in the autocorrelation functions-----!- 3 .. 

Fig. 13 shows the density autocorrelation functions for 
different states on the isochore <f> c — 0.40 approaching 
the attractive glass transition, for different wavevectors: 
at the pre-peak (qa = 1.07) and two higher wavevectors. 
The correlators in Fig. [3] have been time-rescaled to give 
collapse in the final a-decay. Whereas such collapse is 
satisfactory for high wavevectors, at the pre-peak such 
a scaling is not possible: only at the final stages of de- 
cay can the curves can be said to overlap, as shown. For 
the master decay curves, observed at high q where satis- 
factory scaling is accomplished, the decay of the density 
correlators is very stretched, and clear plateaus between 
the short-time /3 relaxation and the terminal a region are 
thus not observed. Nevertheless, the values of the non- 
ergodicity parameters f q , estimated from the heights of 
the ill-defined plateaus that can be seen, are clearly much 
bigger than those of the glass transition in hard spheres—^. 
This was correctly predicted by MCT—. These results are 
qualitatively similar to, but quantitatively different from, 
the incoherent autocorrelators studied previously—-. 

The decay from the plateau of the density autocorre- 
lation function can be described by the von Schweidler 
expression, 

*,(*) = /, - h q (t/r) b [1 - k g (t/r)»] + 0(t 3b ) (5) 

where f q is the nonergodicity parameter, b is the von 
Schweidler exponent, and h q and k q are amplitudes. The 
scaling presented in Fig. [3] shows that these parameters 
are almost the same for all states approaching the transi- 
tion, consistent with their state independence as required 
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FIG. 4: Density correlation functions for <f> c — 0.40 and 
<j> p = 0.42 for qa = 1.07, 3.9, 6.9, 9.9, 15, 20, 25 and 30. 
The fittings using eq. (HJ are shown by the dashed lines. The 
von Scweidler exponent, b, is the same for all q, b = 0.37. 
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FIG. 5: Coherent and incoherent nonergodicty parameters 
(solid line and thick points, respectively). The structure fac- 
tor for 4> p = 0.42 is also included and the Gaussian approxi- 
mation fitting to fg (broken line) . 



by the cv-scaling of MCT, with a change only in the time 
scale r. (This applies for wavevectors that are not too 
low.) A similar expression can be used to describe the 
self part of the intermediate scattering function^, (t) , 
with a different set of parameters /* , h s q and k s q . 

Expression (JSJ is fitted to both the coherent and inco- 
herent density autocorrelation functions; the fittings to 
the coherent correlation functions are presented in Fig. 
0] for various wavevectors. In all the fittings, the von 
Schweidler exponent was b = 0.37, independent of the 
state or the wave vector. This figure shows that the de- 
cay of the density correlators indeed can be described 
using the von Schweidler decay. 

The coherent and incoherent nonergodicity parame- 
ters, f q and /* respectively, and the amplitudes h q and 
hq are presented in Fig. \5\ with the structure factor for 
the state 4> p — 0.42. Whereas the incoherent nonergodic- 
ity parameter f* decays monotonically from 1 to with 
increasing q, the coherent analogue f q oscillates visibly, 
in phase with S q , for qa < 10. At higher qa, both the 
coherent and incoherent nonergodicity parameters decay 
together, without oscillations. These results compare 
nicely with the theoretical predictions 6 -, and show that 
at high q the collective and self dynamics are equivalent. 
But note that for short range attractions it is these high 
q density fluctuations that dominate the arrest transi- 
tion within MCT calculations. This is unlike the case of 
the repulsive glass, where the dominant fluctuations are 
around the main peak in S(q), where the largest differ- 
ences between the incoherent and coherent nonergodicity 
parameters arise. 

The amplitudes h q and h s are shown in the inset in 
Fig. El The time scale r in eq. iJSJ is needed in order 
to get absolute values of these amplitudes, but it cannot 
be accessed from the simulations. We have set r to fulfil 
<$> q (r q ) = fq/e for qa — 9.9. Then h s q presents a maximum 
at qa w 12, whereas h q oscillates at low wavevectors and 
coincides with h q for large q. Contrary to f q , the oscilla- 



tions in the amplitude h q are out of phase with respect to 

S(q), as predicted by MCT for other systems 36 (forming 

repulsive glasses). 

A Gaussian approximation to the van Hove function 

can be used2£ to extract information from ff. In this 

•' i 

simple approximation, for low q, ~ exp{— q 2 (<5r 2 )/6}, 
where (Sr 2 ) is the mean squared displacement. Thus, 
the incoherent nonergodicity parameter yields the lo- 
calization length, n, via f* = exp{— q 2 rf/6} (dashed 
line in Fig. [SJ. The localization length so obtained is 
n = 0.129a (rf = 0.0168a 2 ), about half the size of the 
interaction range, and much smaller than the localization 
length in the repulsion driven glass transition (which is of 
the order of the Lindemann distance, rf ~ 0.1a 2 ). This 
finding shows that the driving mechanism of the arrest 
transition in our system is bond formation, i.e. density 
fluctuations at high q. 

We now define T g , the scaling time to collapse the a- 
decay of the density correlator at a particular wavevector 
q, by the equation $ 9 (r 9 ) = f q /e. The results are plotted 
in Fig. HJ1 According to MCT, all the T q are governed 
by only one q-independent time scale of the a-decay, r 
in expression (J5J, which diverges at the transition point 
cf>p with a power law, while the self diffusion coefficient 
vanishes with a similar form: 

r~(«^-^p A)~(<^-0p) 7 (6) 

The exponent 7 is related to the von Schweidler expo- 
nent b within MCT, leaving only one single parameter de- 
pendent on the interaction details. The transition point 
(ftp has been determined previously using the long-time 
scaling of the incoherent correlation functions^ 3 -. Fig. 
shows that, with <f>p = 0.4265, T q indeed follows power- 
law divergences at all wavevectors, even for qa = 1.07, 
where scaling was not found. The exponents for high 
wavevectors lie in the range 3.36—3.53 without any trend, 
and 7 = 3.08 for qa — 1.07. The time scales of the inco- 
herent correlation functions also diverge following power 
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FIG. 6: Time scale T q at g = 1.07, 3.9, 9.9, 20 and 30 from top 
to bottom, and diffusion coefficient (xlO ), Do, as a function 
of the distance to the transition. The dashed lines represent 
the power-law fittings according to ©. 



laws, with exponents in the range 13 3.1 — 3.25. Using the 
value of the von Schweidler exponent, b — 0.37, the MCT 
relations yield 7 = 3.44, which lies well inside the range 
of observed exponents for the divergence of r q . The sim- 
ilarity between all of the exponents found for the time 
scale shows that, as MCT predicts, the function T q ((f> p ) 
factorizes as T(4> p )uj(q), where u)(q) is a function only of 
q. Closer to the transition point, at (j> p — 0.425, the time 

as 



scale deviates from the MCT power-law behaviour^, 
also observed in other glass forming systems. 

The diffusion coefficient, as obtained from the long 
time slope of the mean squared displacement, is also 
shown in Fig. El A power-law decrease is observed 
for Dq with the same (j>^ , but with a smaller exponent, 
7 = 1.23. Differences between these exponents (7 deter- 
mined from t or Do) have also been observed in simula- 
tions of other systems, such as Lennard-Jones particles^, 
hard spheresS^, polymer melts 3 ^ or silica^, and imply 
the breakdown of the Stokes-Einstein relationship close 
to the glass transition. Alternatively, similar 7 exponents 
can be recovered for Dq and t, so long as different transi- 
tion point values <j>^ are allowed in each fit; however this 
also represents a deviation from MCT predictions. Note 
that in any case, activated processes close to the transi- 
tion will cause deviations from the MCT power-law pre- 
dictions, restoring ergodicity beyond the transition point 
estimated by fitting to the power law region. 

The new results reported here, as well as previous 
analysis of the incoherent correlation functionsiSiiii, show 
that most (though not all) of the predictions drawn from 
the mode coupling theory are confirmed by our simu- 
lations. However, the presence of both structural and 
dynamical inhomogcncitics would lead us to a different 
view of the system, less homogeneous and well behaved, 
which might not be reflected in averaged properties such 
as the coherent and incoherent correlators. 
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FIG. 7: Distribution function P(8r 2 \t) of squared displace- 
ments for the states labeled, equispaced in time t. (These 
times are chosen such that (Sr 2 ) « la 2 for the shortest t and 
(8r 2 ) ~ 50a 2 for the largest; the actual t values therefore vary 
with <f> p .) As time evolves the distribution moves to the right. 
The thick blue lines correspond to (Sr 2 ) » 10a 2 , where the 
distinction between fast and slow particles is made in the fol- 
lowing figures. The dashed vertical line shows the squared 
displacement of a particle in the frozen environment. 
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FIG. 8: Number of fast (open circles) and slow (closed circles) 
particles, and percolation probability (crosses) as a function 
of the polymer fraction. 
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TABLE I: Mean number of neighbors for the whole system, 
n„«, for the fast particles rif and slow particles n s , fast neigh- 
bors of fast particles, n// and slow neighbors of slow particles, 
n ss . Note that the simulations studied here are longer than 
those presented previously--, making the results here more 
reliable. 



B. Dynamical Heterogeneities 

In Fig. [7| the probability distribution P(Sr 2 ;t) of 
squared displacements 5r 2 measured over a time interval 
t is presented (for various t values) , in a series of differ- 
ent states approaching the gel transition. (Curves with 
the same colors in different plots correspond to similar 
mean squared displacements.) In a homogeneous fluid, 
far from any arrest transition, this distribution is single 
peaked at all times; the peak moves to larger distances 
and broadens, with both mechanisms controlled by the 
diffusion coefficienti^. This behavior is seen at the lowest 
polymer fraction studied here (upper panel in Fig. [TJ). 

However, as the polymer fraction is increased, and 
the gel transition approached, the distribution becomes 
wider and, at <p p = 0.42, two peaks become clearly visi- 
ble, indicating two populations of particles with different 
mobilities (at least, on the time scales set by the inter- 
vals t). In fact, the peak corresponding to 'slow' parti- 
cles can already be noticed at (f> p = 0.40 as a shoulder, 
and at <p p = 0.41 as a peak at lower displacements than 
the main peak, but still moving to the right. The slow 
peak at 4> p — 0.42 is, however, at distances of the or- 
der of the localization length, showing a set of particles 
stuck at that distance. Note that the main peak, corre- 
sponding to a population of fast particles, is at squared 
distances three orders of magnitude larger than this by 
the final dataset. Even at short times, when (8r 2 ) w la 2 , 
a shoulder at large distances indicates the existence of a 
population of fast particles^. 

For high polymer fractions, a minimum develops at 
5r 2 ~ a 2 . Once this happens, it is a meaningful approx- 
imation to describe the system with a two population 
model in which fast and slow particles are distinguished. 
At lower polymer fractions, however, this distinction is 
less clear, and only indicative. 

The evolution of the distributions presented in Fig. [7| 
indicates that, as 4> p is raised, not only does the dis- 
tinction between fast and slow particles become clearer, 



but also the fraction of slow particles increases^. This is 
studied in Fig. [SJ where the fraction of slow and fast par- 
ticles is presented as a function of the polymer fraction 
for a fixed time, t* , chosen such that (Sr 2 (t*)) 10a 2 . 
(The distinction between fast and slow particles is made 
by slicing the distribution at 5r 2 — a 2 ; see above.) Inter- 
estingly, at the highest 4> p , the slow particles amount to 
65% of the system, and the square displacement at the 
peak in the distribution describing them increases only 
to ~ 0.1a 2 on the time scales studied; this is in effect, a 
population of particles that remain stuck, very close to 
their initial spatial locations. 

However, the collective motion of the particles, even of 
the slow particles, can be observed by studying the mean 
squared displacements of one single particle moving in 
an environment of frozen particles. (In the simulations, 
only that particle is moved, using the configuration of 
the equilibrated system at the desired stated.) This dis- 
placement is constant up to very long times and is shown 
in the lowest panel of Fig. by a vertical dashed line. 
Its value is much smaller than the peaks of both the fast 
and the slow particles, indicating that although the slow 
particles are practically fixed, their cooperative motion 
allows them to explore longer distances than their cage 
size would be if the environment were frozen. 

We now extend further the analysis of fast and slow 
distributions, aiming to understand better the nature of 
these two populations. In table [I] the mean number of 
neighbors of the fast and slow particles is presented, along 
with the average number of neighbors for the whole sys- 
tem. ('Neighbors' are defined as particles connected by a 
bond, i.e., lying within range of the attractive part of the 
potential.) The slow particles have more neighbors than 
average for all states, whereas the fast particles have less 
neighbors. Thus, we may conclude that the slow parti- 
cles are in the inner parts of a bonding network in which 
both types of particle participate (the 'skeleton' of the 
network), while the fast particles are in the outer lay- 
ers of the network (its 'skin'). As the gel transition is 
approached the difference between the two types of par- 
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Time 

FIG. 9: Fraction of slow particles that were slow in the first 
correlator as a function of time and the correlator index (in- 
set). See text for the latter. From left to right: 4> p = 0.39, 
0.40, 0.405, 0.41, 0.415 and 0.42. The vertical dashed lines 
indicate t*. 



tides is more important, as reflected in their number of 
neighbors, thus implying that this picture of skeleton and 
skin is more appropriate closest to the gel. Note that the 
number of neighbors of the fast particles decreases as 4> p 
increases. Also in table[I]the number of fast neighbors of 
a fast particle, riff, and slow ones of a slow particle, n ss , 
are given for all states. Supporting the simple picture of 
skeleton and skin, at high polymer fractions most of the 
neighbors of a slow particle are slow, and most of those 
of a fast particle are fast. Using the number of neighbors 
of every particle as a measure of the local density, this 
implies spatial correlations for both populations estab- 
lishing further the link between structural heterogeneity 
and dynamic heterogeneity 1 ^. 

Since the slow particles are mostly bonded to each 
other, and they can amount to more than half the system, 
it is possible that the slow particles percolate, forming a 
long lived elastic structure: truly a 'skeleton' of the gel. 
The percolation probability for the slow particles (i.e. the 
probability of finding a cluster of slow particles, percolat- 
ing in all three dimensions, for a given realisation of the 
system of the size that we use) is presented in Fig. [H]as a 
function of the polymer fraction. This probability rises at 
(j)p » 0.40, and percolation is systematically observed for 
<\> v > 0.41. Interestingly, percolation of the slow particles 
sets in over just the same range of <pp as marks the onset 
of a clear distinction between the slow and fast particles 
in the distribution of squared displacements (see Fig. [TJ. 
However, in the following we show that any link between 
these effects is somewhat indirect. 

The picture of skeleton and skin presented above re- 
sembles the channel diffusion of sodium atoms in silica 
melts 41 ^. Channels in the Si-0 network allow diffusion 
of the sodium atoms, resulting in a much higher mobil- 
ity of the latter than for silicon or oxygen. The different 
dynamics is caused by the differences in the interactions, 
but this also results in structural heterogeneities, as in 



the case of gelation. However, in our case, there is no in- 
trinsic difference between the two types of particles, fast 
and slow, and thus, in colloidal gelation, the stability of 
the two populations is not guaranteed, and an exchange 
between them may be ocurring. 

This exchange is addressed in Fig. [§1 where the 'per- 
sistence' of slow particles is presented. We define this 
persistence as follows 14 . We calculate a series of corre- 
lators with different initial times. A new correlator is 
initiated (alongside its predecessors) once the previous 
one has evolved to the point where (Sr 2 ) = a 2 ; this en- 
sures quasi-independence of one correlator from the next 
(because a 2 ^> rf). The fast or slow character of every 
particle is decided after time interval t* for a given cor- 
relator, with t* chosen so that {Sr 2 (t*)) ps 10a 2 . The 
persistence is denned as the fraction of particles, labelled 
as slow in a particular correlator, that remain slow when 
the classification is redone for a subsequent correlator. 
In the inset of Fig. [§1 we show the persistence of the 
slow particles as a function of the correlator index; the 
main panel shows this as a function of elapsed time be- 
tween correlators. (The inset effectively serves as a time- 
rescaling using the mean squared displacement.) As ex- 
pected, the higher the polymer fraction, the slower the 
exchange between the two populations, and the slower 
the decay of the persistence of the slow particles. At low 
polymer fractions, most of the slow particles are reclassi- 
fied as fast on a timescale of a few t* (a few correlators) 
and vice versa; but at high cf> p , only a small fraction of 
slow particles become fast (presumably by escape from 
the percolating bonded network that they form) even at 
times that are an order of magnitude larger than the dy- 
namical relaxation times r q set by the a relaxation time 
r (see Fig. |SJ). These results are quite similar to those 
reported previously for the persistence of fast particles 1 ^. 

The vertical dashed lines in Fig. [§]mark the time when 
the mean squared displacement is equal to 10a 2 , which 
is when the classification into fast and slow particles is 
made. This plot shows that the fraction of slow particles 
that become fast decreases with increasing <f> p , although 
the absolute number of particles is between 90 and 150 
for all the states. In the following we will analyse the 
dynamics of both sets of particles, fast and slow, on the 
basis of a simplified assumption that there is no exchange 
between them. Below t* , this assumption is fairly accu- 
rate, but beyond t* the analysis becomes qualitative only. 

In Fig. we present the coherent intermediate scat- 
tering functions at the wavevector of the prepeak in S(q), 
representative of structure at the length scale of any in- 
cipient microphase separation. The contributions from 
the fast particles (upper panel) and slow particles (lower 
panel), are shown, as well as the globally averaged func- 
tions (inset) 50 . The time t* coincides closely with the 
kink in the curves, signalling the new mechanism enter- 
ing the dynamics, i.e. exchange between the fast and slow 
populations. The contributions from the fast particles 
decay in all cases to below 0.2 before t* is reached, while 
the contributions from the slow particles decay no lower 




FIG. 10: Intermediate scattering function for <f> p = 0.39, 
0.395, 0.40, 0.405, 0.41, 0.415 and 0.42 from left to right. 
Contributions from fast-fast particles (upper panel), slow-slow 
particles (lower panel) and total function (inset). 



than 0.85 by the same time. Beyond t* , the contribution 
from the slow particles decays much faster, possibly due 
to the conversion of some slow particles to fast. Similar 
behaviors were observed in the incoherent functions^, 
showing that the slow particles are tightly bonded to the 
network and their escape occurs only at very long times. 

It should be noted that the contribution from the fast 
particles to both the coherent and incoherent functions 
decays slower for high polymer fractions, showing that 
these particles also feel the proximity of the transition 
(their bonds are longer lived). Thus, as the gel is ap- 
proached, the slowing down of the system is due to all 
particles in the system, and not exclusively to the slow 
ones, whose fraction is increasing also. Also, it is in- 
teresting that no special feature is visible when the slow 
particles percolate, at (j> p = 0.41. This is evidence against 
any direct link between percolation of the slow particles 
and the onset of a clear classification (via the bimodality 
of P(Sr 2 ;t)) into fast and slow populations. 

These features are once again consistent with a skele- 
tal network of stuck particles, surrounded by a skin of 
movable particles. The network is very stable in time, 
and relaxes on a time scale much larger than that of au- 
tocorrelators averaged over the whole system. However, 
the above results do not show whether this relaxation is 
a structural one, or caused by single particles leaving the 
network (thus converting slow particles into fast ones). 
Indeed, our definition of fast and slow particles is based 
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FIG. 11: Environment correlation function for the same 
states as Fig. 1101 Contributions from the fast neighbors of 
fast particles (upper panel), slow neighbors of slow particles 
(lower panel= and total function (inset) are presented. 



only on the squared displacement of every particle, and 
thus it does not distinguish between structural relaxation 
and single particle motions. Therefore, taking the distri- 
bution of squared displacements at different times (dif- 
ferent t*) does provide qualitatively the same results, i.e. 
the correlation function from the slow particles hardly 
decays until t* , when a rapid decay starts (only the frac- 
tion of fast and slow particles varies for different t*). 

To clarify the origin of the long-time relaxation of the 
structure of the slow particles, we have calculated an en- 
vironment correlation function, $ e (i). By calculating the 
fraction of neighbors at given time, t, that were neigh- 
bors at time t = 0, we quantify how the environment 
of the particles changes^. This function is similar to 
the cage correlation function, introduced by Rabani et 
al— or the bond correlation function defined by Luzar 
and Chandler—, but the finite range of the attraction in 
our case provides an unambiguous way to define neigh- 
bors. In Fig. ^2 we present & e {t), and its contributions 
from the fast neighbors of fast particles, and the slow 
neighbors of the slow ones. These contributions are the 
dominant ones close to the gel, according to the distribu- 
tion of neighbors presented in Table [I] The data in Fig. 
1111 extend those reported previously^ in both time and 
composition. 

As expected the total environment correlation func- 
tion decays more slowly at high polymer fractions, sig- 
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FIG. 12: Coherent (continuous lines) and incoherent (broken 
lines) density correlation functions for <j> p — 0.42, and for 
qa = 1.07 (upper panel) and qa = 9.9 (lower panel). From 
top to bottom, contribution from the slow particles, whole 
system, and contribution from the fast particles. 

nailing the slower relaxation of the environment of the 
particles. As already observed in Fig. ^3 the slow parti- 
cles hardly change their environment before t* , whereas 
the fast ones retain less than 20% of their neighbors on 
this time scale; this is again consistent with our 'skele- 
ton and skin' description of the system. At larger times, 
however, the neighborhood of the slow particles finally 
relaxes (see Fig. 9) implying that the particles change 
their neighbors, thus escaping the network. This seems 
to suggest that the relaxation of the network of slow 
particles observed at long times occurs in a dissolution- 
like way, rather than by collective structural relaxation 
(which would allow slow regions to become fast without 
reassignment of neighbors). Such a conclusion must be 
treated cautiously though, since our definition of neigh- 
bors (via bonding) is a strict one: a particle could alter its 
neighbors without structural reorganization on the cage 
scale or above. 

Finally, we study the behaviour of the system at short 
distances by means of the density correlation function at 
high wavevector, qa = 9.9. In Fig. ^|both the coherent 
and incoherent intermediate scattering functions are pre- 
sented for qa = 1.07 (upper panel) and qa = 9.9 (lower 
panel), at <\> v — 0.42. Again the contributions from the 
fast and slow particles to both correlation functions are 
shown. The correlation functions at this higher wavevec- 
tor show that the fast particles are very mobile, but the 



slow ones are more trapped. However, as Fig. [7| already 
showed, they are not completely stuck, and a significant 
decay of Q q is observed before t* for the slow particles 
(note that this contribution at low q decays only to 0.95). 
The bond correlation functions, not shown, indicate that 
40% of the bonds between slow particles break before t* . 
Thus, the decorrelation at high q implies bond breaking, 
and not only structural relaxation below t* , although the 
environment of the particles does not change significantly 
(see Fig. [nj. 



IV. CONCLUSIONS 

In this paper and its predecessorsi2ii2ii4 we have pre- 
sented the analysis of states close to gelation, in a system 
with short range interactions. Recent MCT results indi- 
cate that gelation at high density (as considered here) 
should be viewed as an attraction-driven glass transi- 
tion; our system has been analysed to test these theoret- 
ical predictions. On the other hand, dynamical hetero- 
geneities, although also observed in other glass forming 
systems, are found to be very pronounced in our system, 
and this might be taken as evidence against the MCT 
picture. We have studied these heterogeneities in some 
detail in an attempt to resolve their character. 

The analysis of our system within MCT shows that 
most of the globally averaged properties, as measured by 
coherent and incoherent correlators, are correctly pre- 
dicted. The a-scaling of the density correlation func- 
tion works for wavevectors that are not too low, and the 
system-universality of the von Schweidler exponent b and 
the dynamic exponent 7, and the predicted relationship 
between these two, are confirmed. The driving mecha- 
nism, namely bond formation, causes differences between 
the attraction-driven glass studied here and and caging- 
driven repulsive glasses: the decay from the plateau is 
more stretched, the nonergodicity parameters are higher, 
and the localization length is shorter. All these trends 
are fully predicted by MCT. However, some significant 
differences from MCT predictions are noticed, such as 
the appearence of different exponents in fits for the di- 
vergence of r q and 1 /Do, and the breakdown of the en- 
sealing at low q. (The latter might be specifically related 
to the repulsive barrier used in our potential.) 

By studying the distribution of squared displacements 
at fixed time in the system, it is possible to recognize 
two populations of particles with different mobility, in 
systems close to the gelation transition. The analysis 
of the correlation functions at different wavevectors (as 
well as environmental and bond correlators) shows that 
the slow particles form a quasi-frozen structure, which 
percolates at high enough tp p . (However the threshold is 
discernably lower than </& which, according to the MCT 
fits, marks the onset of gelation itself.) The fast particles 
are in the surface or 'skin' of this cluster, and can escape 
more easily than the slow particles in the inner parts 
of it (the 'skeleton'). There is a fairly clear correlation 
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between the classification into fast and slow particles and 
the presence of structural heterogeneity as quantified in 
the local statistics of bonded neighbors. 

The slow particles, on the other hand, are not com- 
pletely stuck, but can break their bonds and move over 
short distances. Despite this, they form a structure which 
hardly changes on lengthscales above that of a single par- 
ticle. (In other words, broken bonds among slow particles 
are replaced by others before significant rearrangements 
can occur at larger length scales.) Only at very long times 
does the skeleton of slow bonded particles relax signifi- 
cantly its shape; and this seems to be due to a very slow 
exchange between fast and slow particles. This relax- 
ation appears to be dissolution-like, in the sense that sin- 
gle slow particles escape to become fast, while elsewhere 
the reverse process occurs. However, evidence for coop- 
erativity in the motion of fast and, more importantly, 
slow particles, is also found: the squared displacements 
of all particles is larger than that of single particles in 
the frozen environment. 

Within our results, we cannot discern whether all of 
the slow particles will eventually change to fast, causing a 
complete relaxation of the structure as one would expect 
in a fluid phase. Although the curves in Fig. [!|] appear 
to saturate in all cases, the states at low <p p must con- 
tinue the exchange, since these are certainly fluid phases 
(see Fig. [7J|. Moreover, even our highest polymer den- 
sity (j) p = 0.42 is also a fluid phase if the gel transition 
is accurately identified^ from the best von Schweidler fit 
((f>p = 0.4265), in which case the same can be expected 
for all our samples. However, the very presence of persis- 
tent structural and dynamical heterogeneity that endures 
through the entire lifetime of our simulation runs makes 
an accurate identification of the gel transition unreliable 
(to say nothing of issues of equilibration as raised in the 
introduction) . This rather complicated physics arises de- 
spite the apparently 'good behavior' of the globally aver- 
aged correlators - which do seem to behave in line with 
MCT predictions for systems close to arrest (but still 
within the ergodic fluid phase). 

The intriguing relation between the onset of percola- 
tion among the slow particles and the onset of a clear 
bimodality in P(Sr 2 ;t) among fast and slow populations 



of particles merits further study. No direct link can be 
established on the basis of our results, but it may be more 
than coincidence. The formation of a percolating skele- 
tal cluster of slow particles at high enough tp p could be 
instrumental in suppressing diffusive contributions that 
would otherwise arise from collective motions of finite 
clusters, thus sharpening considerably the bimodality of 
the population-level dynamics. However, one might ex- 
pect this process also to imply a visible feature in the 
behavior of the correlation functions at low qa, which is 
not observed. 

The somewhat complicated dynamics exhibited by our 
systems may caution against too sharp a dichotomy be- 
tween different modelling strategies for glasses, at least 
when applied to colloidal suspensions^. If MCT is taken 
as a priori incompatible with dynamic heterogeneity of 
any kind, then our results suggest that MCT's recent 
quantitative prediction of the experimental behavior of 
attracting colloids must in part be fortuitous. And in 
fact, qualitative retrodiction of some of these results 
has recently been achieved using facilitated dynamics 
models 4 -. But note also that very recent MCT work 
claims anyway to encompass DB»£; and a strong case has 
been made for a diverging dynamical lcngthscalc with as- 
sociated fluctuations within MCT—. 

In any case the DH we observe is of discernibly dif- 
ferent character to that seen in repulsion-driven glasses. 
The latter has formed part of the backdrop to recent 
developments in lattice-based models of DH involving fa- 
cilitated dynamics^. A fuller understanding of DH in 
dense colloids with short-range attractions must include 
some of these same elements, but may also involve an 
interplay between arrest and incipient phase separation 
(or, in our simulations, incipient microphase separation) 
which has long been argued to play a strong role in col- 
loidal gelation at low density22^. 
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